Order of the phase transition in models of DNA thermal denaturation 
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We examine the behavior of a model which describes the 
melting of double-stranded DNA chains. The model, with 
displacement-dependent stiffness constants and a Morse on- 
site potential, is analyzed numerically; depending on the stiff- 
ness parameter, it is shown to to have either (i) a second order 
transition with i/± — —/3 = 1, Mi = 7/2 — 2 (characteristic of 
short-range attractive part of the Morse potential) or, (ii) a 
first-order transition with finite melting entropy, discontinu- 
ous fraction of bound pairs, divergent correlation lengths, and 
critical exponents vx = — P = 1/2, V\\ — 7/2 = 1. 

PACS numbers: 05.70.Jk, 63.70.+h, 87.10+e, 05.70.Fh 



Early models of DNA thermal denaturation, i.e. the 
separation of the two strands upon heating Q, were 
based on an Ising-like description in which a base-pair 
was either closed or open; the relative tendencies of the 
system to occupy one of the two states were introduced 
explicitly, in terms of temperature-dependent free en- 
thalpies. As a consequence, although a judicious choice 
of such enthalpies has proved useful in describing some 
aspects of DNA denaturation j|] , understanding of this 
remarkable one-dimensional cooperative phenomenon in 
terms of standard statistical mechanics -i.e. a Hamilto- 
nian model with temperature-independent parameters- 
remained an outstanding problem. 

More recent research has emphasized the role of the 
large amplitude fluctuations that precede the transition 
and the intrinsically nonlinear mechanisms which are 
needed to describe such fluctuations |!|. In such models 
the status of a base pair is characterized by the distance 
between the two bases. An on-site asymmetric poten- 
tial with a flat part at large values of the displacement 
emulates the tendency of the pair to "melt" at high tem- 
peratures, as thermal phonons drive the particles outside 
the well and towards the flat portion of the potential. 

In the original version ("type I" model) coupling be- 
tween successive base pairs is harmonic [Q; the result- 
ing path to the melting instability appears smooth; this 
is at variance with the sharp features of the transition 
observed experimentally. A generalization ||^| of the 
model to include displacement-dependent stiffness con- 
stants, i.e. "stacking parameters" which describe the 
coupling between successive base pairs, has revealed a 
dramatic sharpening of the transition. The predictions 
of the latter, ("type II") model have been compared 



successfully with experimental results 0. Furthermore, 
investigations of heterogeneous DNA have shown that 
the model yields features of multistep melting similar 
to those observed in experiments ||. In fact, the au- 
thors of Ref. H have, in passing, pointed out the formal 
analogies between the melting instability of homogeneous 
DNA and other continuous phase transitions, e.g. wet- 
ting of a one-dimensional interface from a substrate, ad- 
sorption of polymers by a solid surface etc.; in addition, 
they have demonstrated that "type II" models generate 
an " entropic barrier" which is largely responsible for the 
narrowing of the transition. Their analysis of the or- 
der parameter led them however to suggest that, in spite 
of the dramatic narrowing of the transition, the criti- 
cal exponent remains unchanged. As a result, the exact 
character of the homogeneous DNA transition remained 
somewhat elusive. 

In this Letter we report on the scaling behavior of type- 
II model of the denaturation of ideal, homogeneous DNA. 
We show that, for values of the stacking parameter used 
in || the typc-II model exhibits a peculiar type of first 
order transition, with a finite melting entropy, a discon- 
tinuity in the fraction of bound pairs (the usual DNA 
observable), and divergent longitudinal and transverse 
correlation lengths ||; as the value of the stacking pa- 
rameter changes, and the range of the "entropic barrier" 
becomes shorter than, or comparable to the range of the 
Morse potential, the transition changes to second order, 
as in type-I models. 

The Hamiltonian of the model is 
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2m 



W(y n ,y n . 



V(y n ) + Dhay n 



(1) 



where m is the reduced mass of a base pair, y n de- 
notes the stretching of the hydrogen bonds connecting 
the two bases of the n th pair and p n — m(dy n /dt). 
Coupling between successive base pairs is described by 
W{y n ,y n -i) = f + {y n -y n _{f. The 

parameter p can take non-zero values in type-II mod- 
els; The choice of this coupling potential is motivated 
by the observation that the stacking energy is a prop- 
erty of base pairs rather than individual bases. The 
effective coupling constant decreases from (1 + p)K to 
K when either one of the two interacting base pairs is 
open, in qualitative agreement with the known proper- 
ties of base-pair interactions in DNA. The third term in 
(pi) stands for an on-site potential which describes the 
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interaction of the two bases in a pair; the Morse poten- 
tial V(y n ) = D(e~ av " — l) 2 has been chosen because 
it has the correct qualitative shape. Finally, the fourth 
term, which describes the effect of a transverse, exter- 
nal stress h, is in fact a mathematical device useful in 
practical calculations. By letting the dimensionless h 
approach zero from above, it is possible to extract the 
scaling behavior near the transition; at the same time, 
since the partition function is now divergence-free at any 
h > 0, a source of criticism of the model on formal math- 
ematical grounds jl0| is removed. The parameters of the 
model are: D = 0.03 eV, K = 0.06 eV/A 2 , a = 4.5 A" 1 , 
a = 0.35 A -1 , m = 300 a.m.u. and the lattice constant 
I = 4.51. 

The thermodynamic properties of (|l|) can be described 
in terms of the eigenvalues and eigenstates of the 
transfer-integral (TI) equation 

/DO 
dye -K{x,v)/k B T = e -e l/kB T (2) 

-oc 

with a symmetrized kernel K,(x, y) = W(x, y) + ^[V(y) + 
V{x)} + ^Y^(y + x); here T is the temperature and k B the 
Boltzmann constant. Details of the numerical procedure 
used in solving (||) have been given in ||; in the present 
study a Gauss-Legendre quadrature formula has been 
used. In the gradient expansion (continuum) approxima- 
tion, valid in the temperature window D < k B T < K/a 2 
JlH , the integral equation (||) can be well approximated 
by the second-order differential equation 

~ X ri 2 % T + U ( X )<t>i = l i<t>i ( 3 ) 

o g(x)a ax z 

where S = v / 2KD/a 2 /(k B T), g(x) = 1 + pcxp(-2ax), 
De, = a + \k B T l2)\Ti{2ira 2 k B T I K), DU{x) = V{x) + 
(fcflT/2) In g(x), and the limit h — ► has been explicitly 
taken. 

Of particular interest are (i) the lowest eigenvalue eo, 
which, in the thermodynamic limit, is equal to the free 
energy per site /, (ii) the ground state (j>o, which deter- 
mines the order parameter a = (y) = J dyy\4>o\ 2 and its 
fluctuations {{8y) 2 ) = J dy(y — a) 2 \4> \ 2 (= the trans- 
verse correlation length in the language of wetting Q ) , 
and (iii) the next-to-lowest eigenvalue, which controls 
the longitudinal correlation length £|| = lkBT/{€i — eo)- 
Computation of the static structure factor S(q, T) = 
J2 n e-np(-iqnl)(Sy n 5yo) , where Sy n = y n - a, requires 
knowledge of the full spectrum; in terms of the ma- 
trix elements M, = (z|x|0), and the differences Aj = 
(e-i - eo)/fcsr, 

S ( q ,T) = j2W h T hAl (]V (4) 

cosh Aj — cos{ql) 
where the ground state is excluded from the summation. 



For p = (type - I case) the asymptotic properties 
of the denaturation instability can be obtained from the 
solution of the " pseudoSchrodinger" Eq.(||) with g = 1 
Q|5j . In brief: as long as 1 > 5 > S c = 1/2, there is 
a single bound state with io/D = 1 — (1 — S c /S) 2 which 
disappears at S = S c , corresponding to a critical tempera- 
ture k B T c = 2\j2KD i 'a; in the critical regime preceding 
the instability, \t\ < l,t e T/T c — 1, the power laws 
2/£|| oc 1*1^11 , a ex \t\ fi and £j_ cx \t\~ VJ - hold, with = 2 
and v± = —/3 = 1. Furthermore, we have calculated the 
structure factor |l6| in the regime qa <C 1, ^ 1, 

%^) = ^^ F |^||), (5) 

where F(x) = (l/4)x" 2 [l - l/cosh(arcsinh(x)/2)]. This 
implies (i) a zero-field isothermal susceptibility \ — 
lini/ w o a(da/dh)T oc |t| -7 , where 7 = 4, and (ii) critical 
correlations (q!;\\ ^> 1), S(q) oc (qla)~ 2+v with r\ = 0. 

It should be noted that the occurrence of a ther- 
modynamic transition in the one-dimensional model (|l|) 
does not imply a violation of well-known theorems: van- 
Hove's theorem [jl7| does not apply, since the Hamilto- 
nian includes an on-site term; furthermore, since there 
is no symmetry-breaking (or domain-wall-like solitons) 
involved, the standard Landau argument against phase 
transitions in one dimension is also inapplicable. 

Numerical results for a type-II model (p = l,a = 
0.35A~1, a/ a = 0.078), obtained by solving the TI equa- 
tion @ for various values of T and h, are summarized 
in Fig. 1 ]lq| . Below T c , the difference between the two 
lowest eigenvalues, which determines £|| as well as the 
singular part of the free energy, is found empirically to 
satisfy the scaling equation 

ei-e = -fsing = D\t\ <f> [r^/ij > ( 6 ) 

with the limiting forms $(x) = <&q + <&ix + $2X 2 + ■■■(% <C 
1) and $(x) cx x 2 ^(x > 1); it follows that v\\ = 1, 7 = 2 
and — j3 = 1/2. The values of the first two critical ex- 
ponents are in very good agreement with those obtained 
directly from our zero-field results for the longitudinal 
correlation length and the susceptibility, respectively(cf. 
Fig. 2). The order parameter and, to a lesser extent, the 
transverse correlation length (also shown in Fig. 2) re- 
veal significant deviations from pure power-law behavior 
- presumably due to strong transients (cf below); results 
are however roughly consistent with i/± = 1/2. 

It is possible to follow the transition from type-I to 
type-II behavior by continuously varying the stacking pa- 
rameter a, at constant p = 1. We have done this using 
numerically evaluated semiclassical (Bohr-Sommerfcld) 
eigenvalues of Eq. [|. The results, shown in Fig. 3, indi- 
cate that a crossover takes place at a/a « 0.5. Smaller 
values of aj a correspond to a longer -yet finite- range of 
the entropic barrier, compared with that of the Morse 
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potential; one obtains the type-II exponent, v\\ = 1. As 
a I a increases and the range of the entropic barrier be- 
comes shorter than that of the Morse potential, the en- 
tropic barrier becomes irrelevant to critical behavior; one 
obtains the type-I exponent, v\\ =2. 

Equation (^) states that the singular part of the zero- 
field free energy depends linearly on the temperature. In 
other words, the value of the exponent v\\ = 1 implies a 
first-order transition. For our parameter set, the corre- 
sponding melting entropy is AS/kB = [Z?/(fcsT c )]<i>o = 
.75 (cf inset in Fig. 4). 

At this point a comment on the fraction of bound pairs 
is in order; it is this quantity which one measures for 
DNA, using UV absorbance, rather than the order pa- 
rameter. In the type-I model, the probability of finding a 
given base pair at an equilibrium distance smaller than b 
(equal to the fraction of bound pairs, with a proper choice 
of b), is given in terms of the incomplete gamma function, 
i.e. P(y < b,T) = 1 — 7 (2<5 - 1, 2Se- ab )/T(2S - 1), and 
approaches zero continuously as T — * T c , independently 
of the choice of b. This is not the case in type II behavior. 
It can be seen in Fig. 4 that there is a step discontinu- 
ity in the fraction of bound pairs; the exact magnitude 
of the step depends on the choice of b, but the disconti- 
nuity appears to be an intrinsic property of the type-II 
ground state wave function; the subtlety lies in the fact 
that, although there is a long-tail which causes the di- 
vergence of a, a finite weight of |</>o| 2 originates in finite 
displacements. This is to be contrasted with type-I be- 
havior, where more and more weight is shifted to infinity 
as the transition is approached. Therefore, the experi- 
mental detection of the fraction of bound pairs in terms 
of the function P(y < b, T), provides accurate informa- 
tion about the true order of the transition; the order 
parameter itself diverges smoothly at the transition and 
might, due to transients which mask the leading-order 
asymptotics, be less suited for a detailed study, even if it 
were readily available by experimental methods. This is 
a fortunate natural coincidence. 

In summary, we have presented a complete scaling 
analysis of a simple model which has been developed in 
order to describe the melting of "homogeneous" DNA. 
In terms of biophysical applications, our results should 
be complemented by the analysis of Ref. M to account 
for the effects of heterogeneity. We feel nonetheless, that 
there is a direct gain from the analysis presented here: 
the melting of a double-stranded chain - a fairly gen- 
eral problem of biologically motivated statistical physics 
- has been shown to be a true thermodynamic transition 
with different types of critical behavior, depending on the 
details of the interaction parameters. It would be inter- 
esting to speculate on whether other varieties of thermal 
biomolecular denaturation (e.g. protein unfolding) could 
be studied in terms of related models of low-dimensional 
phase transitions. 
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FIG. 1. Type II critical behavior: the dependence of 
(ei — eo)/D\t\ on the scaling variable h/\t\ 3 ^ 2 is shown for 4 
different values of h and a range of temperatures; the dotted 
line marks the slope 2/3. 
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FIG. 2. Zero field results in the type-II case. The longi- 
tudinal correlation length £|| (circles, left axis), the square 
root of the susceptibility \ (triangles, left axis), the order 
parameter a (diamonds, right axis) and its root-mean-square 
fluctuations £± (squares, right axis) as a function of \t\. The 
dashed and dashed-dotted lines correspond to vu =7/2 = 1 
and vx = 1/2, respectively. 
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FIG. 3. The difference between the lowest eigenvalue of Eq. 
|§] and the bottom of continuum band, as a function of the re- 
duced temperature, calculated numerically for various values 
of the ratio a/ a, within the Bohr-Sommerfeld quantization 
scheme. The asymptotic slope, which gives the critical expo- 
nent i/ii, changes from 1 to 2, as the stacking parameter varies 
from type-II (a/a < 1/2) to type-I (a/a > 1/2) values. The 
dashed and dotted lines have slopes of 1 and 2, respectively. 
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FIG. 4. The fraction of bound base pairs P(y < 2A, T) as 
a function of the T/T c for the type I (triangles) and II (dia- 
monds). Inset: the entropy S(T) /Nks) (symbols as in main 
fig.); the length of the double arrow represents the estimate of 
the melting entropy obtained from the scaling Eq. (^). The 
solid lines are "guides to the eye". 
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